#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

/*
 * 拉格朗日插值法：
 * 参数一：插值点数组，第二维下标为0的是X值，下标为1的是对应的Y值
 * 参数二：插值点个数
 * 参数三：需要模拟的点
 */
double LagrangeInterpolation(double TheInterpolationPoints[][2],int NumOfPoint,double SimulationPoint);
/*
 * 多项式相乘/加：
 * 参数一：相乘的第一个多项式的系数
 * 参数二：相乘的第一个多项书的项数
 * 参数三：相乘的第二哥多项式的系数
 * 参数四：相乘的第二个多项式的项数
 * 参数五：最终结果多项式的项数（传出参数）
 * 返回值是最终结果多项式系数数组的首地址
 * 注：多项式数组下标即为该项自变量的级数，值是该项的系数
 */
double *PolynomialMultiplication(const double Polynomial1[],int N1,const double Polynomial2[],int N2,int *FN);
double *PolynomialAdd(const double Polynomial1[],int N1,const double Polynomial2[],int N2,int *FN);

//用于打印多项式
void printPolynomial(double Polynomial[],int N);

int main() {

    //多项式相乘测试
//    double P1[] = {1,2,3};
//    printPolynomial(P1,3);
//    double P2[] = {1,2};
//    printPolynomial(P2,2);
//
//    int FN;
//    double *Final = PolynomialMultiplication(P1,3,P2,2,&FN);
//    printPolynomial(Final,4);
//    printf("%d\n",FN);

    //测试样例一：
//    double TheInterpolationPoints[4][2] ={{-1,-2},{1,0},{3,-6},{4,3}};
//    double  Ans,SimulationPoint = 7;
//    Ans = LagrangeInterpolation(TheInterpolationPoints,2,SimulationPoint);
//    printf("\nSimulationPoint:%lf\tSimulationValue:%lf\n",SimulationPoint,Ans);

    //测试样例二：
//    double TheInterpolationPoints[3][2] ={{2,0.5},{2.5,0.4},{4,0.25}};
//    double  Ans,SimulationPoint = 3;
//    Ans = LagrangeInterpolation(TheInterpolationPoints,3,SimulationPoint);
//    printf("\nSimulationPoint:%lf\tSimulationValue:%lf\n",SimulationPoint,Ans);

    //测试样例三：
//    double TheInterpolationPoints[42][2] ={{4,2},{9,3}};
//    double  Ans,SimulationPoint = 7;
//    Ans = LagrangeInterpolation(TheInterpolationPoints,2,SimulationPoint);
//    printf("\nSimulationPoint:%lf\tSimulationValue:%lf\n",SimulationPoint,Ans);

    //例题：
    double TheInterpolationPoints[6][2] ={{0.4,0.41075},{0.55,0.57815},{0.65,0.69675},{0.80,0.88811},{0.9,1.02652},{1.05,1.25386}};
    double  Ans,SimulationPoint = 0.596;
    Ans = LagrangeInterpolation(TheInterpolationPoints,6,SimulationPoint);
    printf("\nSimulationPoint:%lf\tSimulationValue:%lf\n",SimulationPoint,Ans);
    return 0;
}

double LagrangeInterpolation(double TheInterpolationPoints[][2],int NumOfPoint,double SimulationPoint){

    double **BasisFunctionArray = (double**)malloc(sizeof(double *)*NumOfPoint);  //基函数所在数组
    for(int i = 0;i < NumOfPoint;i++){
        double *temp = (double*)malloc(sizeof(double));
        *temp = 1;  //初始化基函数
        BasisFunctionArray[i] = temp;   //初始化基函数所在数组
    }

    for(int i = 0;i < NumOfPoint;i++){

        //构建基函数
        int NTerms;     //每个基函数多项式相乘后总共的项数
        double DividendKey=1;
        for(int j = 0;j < NumOfPoint;j++){
            if(j != i){     //拉格朗日插值法
                DividendKey/=(TheInterpolationPoints[i][0]-TheInterpolationPoints[j][0]);   //每个基函数最后所除以的因子
                double temp[2] = {-TheInterpolationPoints[j][0],1}; //每个单独的小函数
                BasisFunctionArray[i] = PolynomialMultiplication(BasisFunctionArray[i],1+j,temp,2,&NTerms);     //和每个小函数相乘
            }
        }
        //除以因子得到最后的一个基函数
        DividendKey*=TheInterpolationPoints[i][1];  //每个基函数构造插值多项式时乘的Y值，顺便就乘在了因子上
        for (int j = 0; j < NTerms; ++j)
            BasisFunctionArray[i][j]*=DividendKey;

        printf("\nThe basis function divided by the factor:\n");
        printPolynomial(BasisFunctionArray[i],NTerms);
    }

    double *InterpolationPolynomial = malloc(sizeof(double )*NumOfPoint);
    memset(InterpolationPolynomial,0, sizeof(double )*NumOfPoint);

    int FN; //最终插值多项式的项数
    for (int i = 0; i < NumOfPoint; ++i) {  //构造插值多项式
        InterpolationPolynomial = PolynomialAdd(InterpolationPolynomial,NumOfPoint,BasisFunctionArray[i],NumOfPoint,&FN);
    }

    printf("\nLagrange interpolation polynomials:\n");
    printPolynomial(InterpolationPolynomial,NumOfPoint);

    double Ans=0;
    for (int i = 0; i < NumOfPoint; ++i) {
        Ans+=InterpolationPolynomial[i]* pow(SimulationPoint,i);
    }
    return Ans;

}

double *PolynomialAdd(const double Polynomial1[],int N1,const double Polynomial2[],int N2,int *FN){
    *FN = N1>N2?N1:N2;  //获取两个多项式中最常的多项式的项数
    //为结果多项式开辟空间
    double *FinalPolynomial = (double *)malloc(sizeof(double )*(*FN));
    memset(FinalPolynomial,0, sizeof(double )*(*FN));

    //双指针，分别指向两个多项式的每一项
    int p1=0,p2=0;
    while (p1<N1||p2<N2){   //只要两个多项式都没有加到头
        if(p1<N1){  //多项式对应次幂的系数相加
            FinalPolynomial[p1] += Polynomial1[p1];
            p1++;
        }
        if(p2<N2){  //多项式对应次幂的系数相加
            FinalPolynomial[p2] += Polynomial2[p2];
            p2++;
        }
    }
    //返回答案多项式
    return FinalPolynomial;
}

double *PolynomialMultiplication(const double Polynomial1[],int N1,const double Polynomial2[],int N2,int *FN){
    *(FN) = (N1-1)+(N2-1)+1;    //获取多项式相乘后，应该有的项数
    //为结果多项式开辟空间
    double *FinalPolynomial = (double *)malloc(sizeof(double )*(*FN));
    memset(FinalPolynomial,0, sizeof(double )*(*FN));

    //多项式每一项逐一相乘，下标就是对应次幂，值就是对应次幂的系数
    for (int i = 0; i < N1; ++i) {
        for (int j = 0; j < N2; ++j) {
            FinalPolynomial[i+j] += Polynomial1[i]*Polynomial2[j];
        }
    }
    //返回结果多项式
    return FinalPolynomial;
}

//用于打印多项式，检测结果
void printPolynomial(double Polynomial[],int N){
    for (int i = 0; i < N; ++i) {
        printf("(%.5lf)x^%d\t",Polynomial[i],i);
    }
    printf("\n");
}